Improved modelling of liquid GeSe 2 : the impact of the exchange-correlation functional 
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The structural properties of liquid GeSe2 are studied by using first-principles molecular dynamics 
in conjuncton with the Becke, Lee, Yang and Parr (BLYP) generalized gradient approximation for 
the exchange and correlation energy. The results on partial pair correlation functions, coordina- 
tion numbers, bond angle distributions and partial structure factors are compared with available 
experimental data and with previous first-principle molecular dynamics results obtained within the 
Perdew and Wang (PW) generalized gradient approximation for the exchange and correlation en- 
ergy. We found that the BLYP approach substantially improves upon the PW one in the case of the 
short-range properties. In particular, the Ge— Ge pair correlation function takes a more structured 
profile that includes a marked first peak due to homopolar bonds, a first maximum exhibiting a clear 
shoulder and a deep minimum, all these features being absent in the previous PW results. Overall, 
the amount of tetrahedral order is significantly increased, in spite of a larger number of Ge— Ge ho- 
mopolar connections. Due to the smaller number of miscoordinations, diffusion coefficients obtained 
by the present BLYP calculation are smaller by at least one order of magnitude than in the PW 
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I. INTRODUCTION 

Disordered Ge n Sei_„ materials feature a large va- 
riety of bonding behaviors as a function of the 
composition^^^ In the liquid state, a delicate inter- 
play between the covalent and the ionic characters sets 
in for increasing values of n. At n=0.33, this results in 
a network system (GeSe2) made of predominant GeSe4 
tetrahedra coexisting with homopolar bonds and defec- 
tive Ge— Se coordinations^ The existence of chemical 
disorder in an otherwise prevailing tetrahedral network 
has been firmly established through the measurement 
of the full set of partial structure factors and pair dis- 
tribution functions of liquid GeSe2 (Z-GeSe2 hereafter). 7 
Early molecular dynamics models based on interatomic 
potentials were unable to predict miscoordinations and 
homopolar bonds£ Very recently, a refined model poten- 
tial for disordered GeSe2 was made available This 
potential reproduces qualitatively the experimental data 
on glassy GeSe2, including the presence of Se— Se ho- 
mopolar bonds. However, Ge— Ge homopolar bonds 
were absent within this description. This proves that 
the explicit account of the electronic structrure in the 
expression of the interatomic forces is crucial to de- 
scribe properly disordered GeSe2 networks. Along these 
lines, two approaches based on density functional theory 
(DFT) have been employed to study Z-GeSe2 via molec- 
ular dynamics. D. Drabold and coworkers have adopted 
a DFT framework based on a nonself-consistent elec- 
tronic structure scheme, the local density approximation 
of DFT and a minimal basis set i n i 12 As an alternative, 
the fully self-consistent evolution of the electronic struc- 
ture described within DFT (i.e. first-principles molec- 



ular dynamics, FPMD in what follows) has been pur- 
sued in the case of /-GeSe2 with the use of plane waves 
and pseudopotentialsi 13 i 14 ' 15 A comparison of the struc- 
tural properties obtained within these two approaches is 
provided in a recent paper for the case of amorphous 
GeSe 2 ^ 

Turning our attention to the FPMD approach and 
to the case of ^-GeSe2, it is worth recalling the indi- 
cations collected through the use of the local density 
approximation (LDA) within DFT. This had the ef- 
fect of producing an atomic structure affected by an 
excessive amount of chemical disorder and homopolar 
bondsJ^ The absence of the FSDP (first sharp diffrac- 
tion peak) in the total neutron structure factor could be 
correlated to the lack of a predominant structural unit 
(the GeSe4 tetrahedron), with comparable percentages 
of Ge atoms two-fold, three-fold, four-fold and five-fold 
coordinated,-- Interestingly, interatomic potentials fea- 
turing formal charges on the Ge and the Se atoms (+4 
and -2, respectively) are able to provide a GeSe2 network 
based on undcfective GeSe4 tetrahedra, together with a 
FSDP in the total neutron structure factor^ In the search 
of a DFT scheme able to recover a network structure fea- 
turing the predominant presence of GeSe4 tetrahedra, the 
generalized gradient approximation (GGA) for the ex- 
change and correlation energy proposed by Perdew and 
Wang (PW hereafter) was adopted^ This choice yields 
a very good agreement with experiments for the total 
neutron structure factor over the entire range of momen- 
tum transfer^ The improvements brought about by the 
GGA in the PW form were found to be due to a better 
account of the ionic character of bonding, as shown in 
Ref. [U through an analysis of the contour plots for the 
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valence charge densities. The larger ionicity of bonding 
introduced by the PW approach manifests itself through 
a larger depletion of the valence charge at the Ge sites 
and a larger accumulation around the Se atoms, the co- 
valency remaining essentially equivalent in the LDA and 
in the PW scheme.— 

Despite this success, detailed comparison of the partial 
correlations revealed the existence of residual differences 
between theory and expriment. The most important of 
these differences concerns the first sharp diffraction peak 
in the concentration-concentration structure factor that 
appears in the experiment but it is absent in our level of 
theory^ In Ref. |l5) these shortcomings were attributed to 
an insufficiently accurate description of Ge— Ge correla- 
tions. This was confirmed by the shape of the calculated 
Ge— Ge correlation function, much less structured than 
its experimental counterpart and by the excessively long 
(15 % more than the experimental value) first-neighbors 
Ge— Ge distances. Longer interatomic Ge— Ge distances 
and less structured Ge— Ge pair correlation functions 
were correlated to an overestimate of the metallic char- 
acter in liquid GeSe 2 . This observation can be taken as a 
guideline for the choice of alternative GGA recipes capa- 
ble of further improving the performances of DFT and, 
in particular, the distribution of the valence charge den- 
sities along the bonds. In what follows, we are interested 
in GGA functionals enchancing a localized distribution 
of the valence electrons at the expenses of a delocalized 
one, intrinsic in schemes inspired by the uniform electron 
gas model, as the PW one. 

With this purpose in mind, we present a new set of 
structural data for short and intermediate range proper- 
ties in Z-GeSe2, based on the GGA scheme after Becke 
(B) for the exchange energy) and Lee, Yang and Parr 
(LYP) for the correlation energ y 18 i 19 This recipe makes 
no assumption on the uniform electron gas character of 
the correlation energy. Therefore, it is a good candidate 
to correct part of the drawbacks of the PW scheme, as 
those arising from an overestimate of the metallic char- 
acter of bonding. The comparative analysis carried out 
with available PW data reveals that short-range prop- 
erties and the diffusion behavior compare significantly 
better with experiments within the the BLYP approach. 
The improvement is less significant for the intermediate 
range properties. 

This paper is organized as follows. In Sec. II, we de- 
scribe our theoretical model. Our results are collected in 
two sections, devoted to real space properties (Sec. Ill) 
and reciprocal space properties (Sec. IV). A summary 
of general considerations on the modelling of disordered 
network-forming material can be found in Sec. V. Con- 
clusive remarks are collected in in Sec. VI. 



II. THEORETICAL MODEL 

Our simulations were performed at constant volume 
on a system consisting of 120 atoms (40 Ge and 80 Se). 



We used a periodically repeated cubic cell of size 15.7 A, 
corresponding to the experimental density of the liquid 
at T=1050 K. We refer to our previous results on l- 
GeSe2 for an extended rationale on the choice of our 
system size. 15 ' 22 The electronic structure was described 
within density functional theory (DFT) and evolved self- 
consistently during the motion^ 

In a series of previous papers, we had adopted the PW 
scheme due to Perdew and Wang as a generalized gradi- 
ent approximationj 13 i 14 ' 15 i 16 i 22 i 35 ' 36 This amounts to go 
beyond the local density approximation (LDA) by us- 
ing an analytic representation of the correlation energy 
e c {p) for a uniform electron gas. This representation al- 
lows for variations of e c (p) as a function of p and the spin 
polarization.— As mentioned in the Introduction, the use 
of the PW scheme improves substantially upon LDA the 
description of both short and intermediated range or- 
der in Z-GeSe2- However, residual deficiences are found 
within the PW when focusing on the Ge— Ge correla- 
tions, in particular at the level of the Ge— Ge interatomic 
distance and pair correlation function, These signatures 
point out an overestimate of the metallic character of 
bonding, implying that the ionicity at the very origin 
of the formation of the tetrahedral GeSe4 coordination 
is undermined by an unphysical accumulation of valence 
electron along the bonds. 

In the search of a GGA recipe correcting these de- 
fects, we resort in this work to the generalized gradi- 
ent approximation after Becke (B) for the exchange en- 
ergy and Lee, Yang and Parr (LYP) for the correlation 
energy. 18 ' 19 It has to be reminded that the record of re- 
liability of the BLYP approach has been firmly assessed, 
as shown in a detailed comparative study of the perfor- 
mances of a large variety of density functional methods. 20 
In Ref. |20|, it was concluded that the BLYP method is the 
DFT method with the best overall performances for prop- 
erties encompassing equilibrium geometries, vibrational 
frequencies and atomization energies of nanostructured 
systems. Furthermore, our choice is motivated by the 
consideration that no explicit reference to the uniform 
electron gas is made in the derivation of the LYP corre- 
lation energy. In the original LYP derivation, a correla- 
tion energy formula due to Colle and Salvetti is recast in 
terms of the electron density and of a suitable Hartree- 
Fock density matrix, providing a correlation energy and 
a correlation potential. 19 i 21 This scheme is expected to 
enhance the localized behavior of the electron density at 
the expenses of the electronic derealization effects that 
favor the metallic character. These effects are built in 
GGA recipes having the uniform electron gas as refer- 
ence system, as the PW one. For this rationale to be 
generally applicable, it would be desirable to understand 
which systems and bonding situations are likely to be 
better described by BLYP than by PW. By focusing on 
multicomponent systems A„B(i_ n ) of concentrations n, 
the BLYP scheme is expected to be more suited than 
the PW one to treat bonding situations characterized by 
a moderate difference of electronegativity (termed A e ; 
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hereafter) among the system components. This is espe- 
cially true for compositions close to the stoichiometry, to 
be intended in this context as the composition at which 
optimal coordination between the species A and B oc- 
curs (for instance GeSe2 within the Ge n Se(i_ n ) family). 
In the case of a large A e ; , the ionic contribution to bond- 
ing is sufficiently large to ensure effective charge trans- 
fer in a way essentially independent on the details of the 
exchange-correlation functional. This is the case of disor- 
dered SiC>2 (A e z=1.54), well described as a corner-sharing 
network within LDAj21 On lowering A e ;, the amount of 
the valence charge density along the bonds becomes non- 
negligible and the relative weight of the ionic and co- 
valent character more delicate to quantify. The case of 
Z-GeSe2 is a prototype of this situtation, being character- 
ized by A e ;=0.54. By choosing BLYP to improve upon 
PW (and LDA) one attempts to mimimize electronic de- 
localization effects that are undesirable in any bonding 
situation characterized by competing ionic and covalent 
contributions. These ideas withstanding, we stress that 
the use of a GGA recipe localizing the valence electron 
density on the atomic sites cannot be legitimated by the 
simple consideration of the A e ; value. Indeed, it depends 
also crucially on the variation of bonding with thermody- 
namic parameters as pressure and temperature. For in- 
stance, in the case of ^-GcSc2 at high temperatures, the 
delocalized character of bonding grows at the expenses 
of the localized one (due to gap closing effects) , thereby 
making much less crucial the use of any functional local- 
izing the electronic valence on the atomic sites£ 

In our work, valence electrons were treated explicitly, 
in conjunction with normconserving pseudopotentials of 
the Trouiller-Martins type to account for core-valence 
interactions j 2 -^ The wave functions were expanded at the 
r point of the supercell on a plane wave basis set defined 
by an energy cutoff of 20 Ry. We found bond lengths 
(do = 3.99 bohr) and vibrational frequencies (u> = 398 
cm -1 ), reproducing the experimental data^ and the cor- 
responding PW results quoted in Ref. [l5| to within at 
most 2 % and 4 %, respectively. 

One configuration extracted from the fully equilibrated 
trajectories obtained in Refs. [HI for liquid GeSe2 with 
E c = 20 was taken as initial set of coordinates for the 
present, new set of calculations carried out within the 
BLYP scheme. We used a fictitious electron mass of 600 
a.u. (i.e. in units of m e a^, where m e is the electronic 
mass and ao is the Bohr radius) and a time step of St = 
0.1 fs to integrate the equations of motion. Temperature 
control is implemented for both ionic and electronic de- 
grees of freedom by using Nose-Hoover thermostatSi 27 i 28 
We carried out simulations at T=(1050±10) K over a 
time periods of 40 ps, with statistical averages taken af- 
ter discarding an intial segment of 2 ps. Since the use 
of the BLYP scheme is aimed at weakening the metal- 
lic character of bonding when compared to the PW one, 
it is worthwhile to check whether the comparison of the 
electronic densities of states for the PW and BLYP case 
sunstantiate this choice. As shown in Fig. [TJ our selec- 



tion of the exchange-correlation functional is legitimated 
by the comparison between the corresponding time av- 
erages (taken at the same temperature T=(1050±10) K) 
of the electronic densities of states. Indeed, the BLYP 
approach is seen to provide a deeper pseudo-gap around 
the Fermi level than the PW one. 



III. REAL SPACE PROPERTIES 
A. Pair correlation functions 

Partial pair correlation functions g a f3{f) are shown in 
Fig. [2j Peak positions and number of neighbors within 
given integration ranges are displayed in Table [T] Among 
the three pair correlation functions, gGeGe(f) is the one 
most affected by the choice of the exchange-correlation 
functional. The BLYP scheme improves upon the PW 
one by yielding a clear first maximum, due to homopo- 
lar Ge— Ge bonds, and a very pronounced first minimum, 
closely reproducing the trends observed in g^Gei 7 ")- ^ n 
9§eGe'( r ) the position of the first peak approaches the 
experimental value (r = 2.45 A, BLYP; r = 2.70 A, 
PW, Ref. M ; r = 2.33 A, Ref. i). Equally favor- 
able is the BLYP prediction of the number of Ge in the 
first-neighbor shell (0.22, see Table [l] to be compared 
with 0.25, Ref.i) The shape of gceG e F ( r ) reproduces the 
shoulder in the main peak occurring at r~3.lA, indica- 
tive of edge sharing (ES) connections among tetrahedra?^ 
The Ge— Se pair correlation function g e Q^g e {T) is charac- 
terized by a prominent main peak and a deep minimum. 
The position of the main peak is slightly displaced toward 
shorter distance, by 0.05 A. The BLYP scheme is able to 
reproduce accurately the height of the first maximum and 
the abrupt decay from the first, sharp maximum down 
to vanishing values. The first shell of coordination has a 
number of neighbors (3.55) in very good agreement with 
experiments (3.50, see Table [I]). Little improvement is 
found for larger distances, both gQeSe P ( r ) anc ^ 9GeSe( r ) 
lacking of the second, small maximum visible in g^s^r). 
Both gseSe P ( r ) an d g P Sse( r ) follows closely the experi- 
mental ffses e (0 for r > 3 A. The present set of calcu- 
lations provides Se— Se correlations very similar to those 
obtained in Ref. [HI, as shown by the close numbers for 
the first coordination shell neighbors (see Table HJ. Ho- 
mopolar Se— Se bonds are found at a distance only 2% 
larger than in the PW case. 

We obtain partial (riG e , ns e ) and average (n) coor- 
dination numbers from the first-neighbor coordination 
numbers, n GeGe , n GeSe , and n Se Se, given in Table [TT1 
The resulting theoretical values are compared to exper- 
imental data in Table [TTJ As pointed out in Ref. |l5j, a 
compensation occurred in the PW case between the un- 
derestimated value of nceGe and the overestimated value 
of n Ge se, leading to a good agreement for n Ge . In the 
present case, the small difference between nce(BLYP) 
and its experimental counterpart is the result of close 
values for each single contributions, i.e. nQ eGe , riGeSe, 
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and nseSe- As a result, the calculated and experimental 
average coordination numbers n differ by only 3.5%. 

Overall, the present calculations improve the short- 
range structure of Z-GeSe2, featuring more structured 
9GeGe( r ) ano - 9 < GeSei r )- ^ n particular, the BLYP ap- 
proach provides much shorter Ge— Ge distances and a 
well defined first shell of Ge neighbors, bringing pair cor- 
relation functions in better agreement with experiments. 



B. Coordination numbers and bond angle 
distribution 

Further insight into the network topology and its sensi- 
tivity to the specific exchange-correlation functional can 
be obtained through n a (l). We define this quantity as the 
average number of atoms of species a /-fold coordinated 
(see Table UTTf . where a are Ge or Se atoms. Consistently 
with the choice made in Ref.Qjl, we here used a cutoff dis- 
tance of 3 A, which corresponds to the first minimum in 
the Ge— Se pair correlation function and well describes 
the first shell of neighbors also for Ge— Ge and Se— Se 
correlations. As a first observation, one notices that the 
BLYP scheme is characterized by a higher proportion of 
Ge-GeSe 3 connections (as much as 23%), contributing 
to a percentage of Ge fourfold coordinated atoms mod- 
erately larger than in the PW case (66% against 61%). 
The increase of Ge— Ge homopolar bonds and Ge-GeSe3 
connections takes place at the expenses of the undefective 
GeSe 4 tetrahedra, lowering from 53.8% (PW) to 41.8% 
(BLYP). These data are a further manifestation of an 
increased number of Ge— Ge first-shell neighbors. The 
distribution of miscoordinated Ge atoms is different in 
the two situations: PW and BLYP favor threefold (/=3) 
and fivefold (1=5) connections, respectively. In particu- 
lar, in Ge-GeSe4 connections a homopolar Ge— Ge bond 
is found to coexist with an adjacent tetrahedral arrange- 
ment. Coordination of Se atoms reflects the increase of 
chemical order occuring within the BLYP scheme. The 
number of twofold cooordinated Se atoms becomes more 
than 10% larger, mostly due to the predominant Se-Ge2 
configuration. A corresponding decrease in the number 
of miscordinations (in particular Se-SeGe2 , lowering from 
13% to 6.4% is noticeable in Table Ell 

In Fig. [3] we show the distribution of the Se— Ge— Se 
(#SoGcSo) and Ge-Se-Ge (6> GoScGc ) bond angles. These 
distributions have been calculated by including neigh- 
bors separated by less than 3 A. Two features are worth 
pointing out. First, the Se— Ge— Se bond angle distribu- 
tion becomes highly symmetric around 109° in the BLYP 
description and has also a higher intensity. This stands 
for an improved tetrahedral order when compared to the 
results of Ref. UM, where the maximum occurred at 103° 
and the distribution has a larger width. Second, two 
distinct peaks are visible in the BLYP Ge— Se— Ge bond 
angle distribution at about 80° and 100° at the place of 
a flat maximum in between the same values. In Ref. [TH 
we showed that the Ge— Se— Ge bond angle distribution 



can be decomposed in two contributions, one associated 
to edge-sharing tetrahedra (for values close to 80°) and 
the other due to corner-sharing tetrahedra (for values 
close to 100°). Similarly, the peaks visible in the shape 
of BLYP Ge— Se— Ge bond angle distribution have to be 
ascribed to these two different connections. Their promi- 
nent appearance stems from a improved tetrehedral or- 
ganization, in line with the analysis of the coordination 
environement for Ge and Se. 

Finally, we compare the number of Ge atoms that 
belong to zero, one, and two fourfold rings. We used 
the counting algorithm based on the shortest-path cri- 
terion first proposed by King and then improved by 
Franzblau i 29 ' 30 ^ 31 Cutoff radii were taken equal to 3.0 A 
for Ge— Ge, Ge— Se and Se— Se interactions. By using 
the ring statistics results, Ge atoms can be termed Ge(0) 
(Ge atoms not belonging to any fourfold ring), Ge(l) (Ge 
atoms belonging to one fourfold ring), and Ge(2) (Ge 
atoms belonging to two fourfold rings). Ge(l) and Ge(2) 
form edge-sharing connections, while Ge(0) encompasses 
not only those Ge atoms involved in corner-sharing con- 
nections (iVGe(CS)) but also some of the Ge atoms form- 
ing homopolar bonds, A Gc _ Gc . In Ref. [TH it was found 
that 55% of the Ge atoms do not belong to fourfold rings 
(Ge(0)), 36% belong to a single fourfold ring (Ge(l)), 
and 9% belong to two fourfold rings (Ge(2)). According 
to the present results obtained with the BLYP scheme 
61% of the Ge atoms are Ge(0), 34% are Ge(l) and 
5% are Ge(2). Therefore, Ge atoms involved in edge- 
sharing connections (iV Gc (ES)), according to PW and 
BLYP calculations, turn out to be Y Go (ES,PW)= 45% 
and Y Go (ES,BLYP) 39%, respectively. 

To obtain 7V Ge (CS), we adopted the proposal of Ref.@, 
i.e. A Gc (CS)= l-Noe(ES)- A Gc „ Go , which holds in the 
absence of extended chains ? By using the results of Ta- 
bled (Ar Gc _ Go (BLYP)=22%, A Go _ Gc (PW)=4%) one ob- 
taines that A r Ge (CS, PW)= 51% and N Ge (CS, BLYP)= 
39%, leading to W Gc (CS, PW)/A Go (ES,PW) = 1.13 and 
Y Go (CS, BLYP)/iV Go (ES,BLYP)= 1 (these values are 
also collected in Table llllj) . The full partial structure 
factor analysis of liquid GeSe2 points toward compara- 
ble numbers for the edge-sharing and corner-sharing sites 
(Y Go (CS, BLYP)/A Go (ES,BLYP)~ 1)1 This is consis- 
tent with the experimental prediction of Ref. 0, where 
the number of edge-sharing and corner-sharing sites were 
found very close. Therefore, the BLYP approach proves 
more suitable to yield the correct relative proportion of 
corner-sharing and edge-sharing tetrahedra. 



IV. RECIPROCAL SPACE PROPERTIES 

A. Faber-Ziman and Bhatia-Thornton partial 
structure factors 

In Fig. [4j we display the comparison between the two 
sets of calculated Faber-Ziman (FZ)22 partial structure 
factors (PW scheme, Ref. UH, and BLYP scheme, present 
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results) and their experimental counterpart*^ In Ref. 
HEl it was established that the performance of PW cal- 
culations were very satisfactory for k values characteris- 
tic of short range properties (k > 2 A -1 ). However, a 
major disagreement existed in the region of the FSDP 
(first sharp diffraction peak), at k <~1. A -1 , in particular 
for the So^ e (k) structure factor and, to a lesser extent, 
for the S§^ e (k) structure factor. Moreover, S^ Ge (k) 
was found to be less structured and slightly shifted to- 
wards smaller wavevectors with respect to the experimen- 
tal curve. For k > 1.5 A -1 a clear improvement is notice- 
able in the shape of SceGe{k)) around the second max- 
imum, that superposes to the experimental result, the 
peak position being located at 2.2 A -1 (PW calculations, 
1.88 A" 1 , experimental value, 2.2 A -1 .-^) The position 
of the FSDP in Sg e a e {k)) is shifted leftward improving 
upon the result of Ref. [jjl in both location and height. 
However, a second spurious peak in the FSDP region 
shows up at larger k values, suggesting that our statis- 
tical accuracy might be further improved by longer runs 
and/or larger simulation cells. In the case of SoeSeik), 
the FSDP height is reduced and the shape of the first 
minimum follows closely the experimental profile, bring- 
ing theory in better agreement with the experimental re- 
sult. As to S^eSe(k), both the present approach and the 
one of Ref. [l5| performs very well over the entire range 
of k values. These pieces of evidence confirm that BLYP 
concur to improve the short-range behavior pertaining to 
Ge atoms, while the changes induces in the intermediate 
range behavior are smaller. 

In this context, it is worthwhile to analyze the com- 
parison between theory and experiments by considering 
the Bhatia- Thornton 3 - 3 , partial structure factors SNN(k) 
(number- number) , S^cik) (number-concentration) and 
Scc{k) (concentration-concentration) (see Fig. [5]) These 
can be obtained by linear combinations of the FZ struc- 
ture factors^ In terms of the Bhatia-Thornton structure 
factors, the total neutron structure factor Sx{k) reads: 

S T (k) = S NN (k) + A[Scc(k)/c Ge c Se -l} 

+BS NC (k), (1) 

where A = c Ge c Se Ab 2 / (b) 2 , B = 2A6/(6), Ab = b Ge - 
bse, (b) = ccebae+csebse, c a and b a denoting the atomic 
fraction and the coherent scattering length of the chem- 
ical species a (6<3e=8.185 fm, &s e =7.97 fm) 7 . This leads 
to coefficients A and B equal to 1.6 xlCP 4 and 0.053, 
respectively. As detailed in Ref. SjvaK^) is a very 
good approximation for the total neutron structure fac- 
tor S T (k), i.e. \S T (k)-S NN (k)\ < 0.015. 

As shown in Fig. both S%%(k) and S§j? p (k) 
are in excellent agreement with experiments, following 
closely the experimental data over the entire range of k 
values. BLYP improves upon PW in the position of the 
FSDP, that is closer to the measured value (0.98 AA _1 , 
Ref. 0, 1.01 AA-\ present results, 1.13 AA" 1 , Ref. 
Il5h and better reproduces the experimental data in the 
range 2 A^ 1 < k < 4 A^ 1 . Similar performances are 
also recorded for the case of <SVc(^)- Most interesting 



is the case of the concentration-concentration structure 
factor S^(k) in the region k < 2. A -1 . In recent 
years, the occurrence of the FSDP in the concentration- 
concentration structure factor S^(k) has stimulated 
intense experimental and theoretical work, in the search 
of its microscopic origins j 34 ' 35 i 36 ' 37 With this purpose in 
mind, the intensities of the FSDP in the Scc(k) have 
been compared for a series of liquid and glasses^ It was 
showed that the FSDP in Scc(k) occurs for moderate 
departures from chemical order, but vanishes cither 
when the chemical order is essentially perfect or for high 
levels of structural disorder. This is exactly the case 
of liquid GeSe2 modeled within the PW scheme, for 
which no FSDP appear in the Scc(k)' 13 ^ 15 Given these 
premises, the question arises on whether the higher 
tctrahedral order exhibited by the present model of 
liquid GeSe2 scheme does result in any improvement 
in the behavior of the FSDP in Scc{k). Although a 
small shoulder has appeared at the FSDP location in 
Sqq YF (fc), the prominant experimental peak remains 
largely underestimated. The results on the partial 
structure factor support the notion that the use of the 
BLYP generalized gradient approximation has a much 
larger effect on to the short-range properties than on the 
intermediate-range ones. 



V. DYNAMICAL PROPERTIES 

Diffusion coefficients are sensitive probes of the struc- 
tural organization in disordered network-forming materi- 
als. On the one hand, in a chemically ordered network, 
the tetrahedra are the main constitutive units, their sta- 
bility strongly reducing atomic mobility. On the other 
hand, departures from chemical order (homopolar bonds, 
miscoordinations) favor atomic mobility since the net- 
work can seek the energetically most favorable arrange- 
ments through bond coordination changes. In a recent 
paper, we have provided a revealing example of the cor- 
relation existing between the values of the diffusion coef- 
ficients and the network structured In Ref. [H we have 
compared structural and dynamical properties of liquid 
GeSe2 as obtained from FPMD and an effective poten- 
tials based the polarizable ionic model (PIM). While a 
sizeable departure from chemical order is found in the 
FPMD model, the PIM structure is highly chemically 
ordered and the GeSe4 tetrahedron is largely predomi- 
nanant among the structural units. As a consequence 
of these drastic structural differences, the diffusion co- 
efficients pertaining to the PIM model are as low as 
- 1 x 10~ 6 cm 2 s _1 at T=3000 K, while they are at least 
ten times larger in the FPMD case at T=1050 K (see 
below). Having established that the tetrahedral order is 
larger in the BLYP model than in the PW one for l-GeSe2 
at T=1050 K, it is of interest to see to what extent this 
is reflected in the values of the corresponding diffusion 
coefficients. 
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The comparison between the calculated statistical av- 
erage of the mean square displacement 



<r 2 {t) >=_(^|r m (i)-r M (0)| 2 } 



(2) 



for both species a, Ge and Se, obtained within the PW 
(Ref.) and the BLYP schemes, is shown in Fig. [SJ In Eq. 
(HI, Tia(t) is the coordinate of the ith particle at time t 
and N a is the number of particles of the species a. The 
diffusive regime is associated with a slope equal to 1 in the 
infinite time limit linear behavior of log <r 2 (t)> vs log t. 
Provided this condition holds, the diffusion constant can 
be obtained as 



D 



< r 2 {t) > 
6t 



(3) 



In Ref. yja, asymptotic values of <r 2 (t)> are 
attained after 5—6 ps to give diffusion co- 
efficients of D Ge =(2.2±0.2)xlCT 5 cm 2 /s and 
D Se =(2.2±0.2)xlO~ 5 cm 2 /s. In the present case, 
establishment of a diffusive regime requires longer time 
intervals and asymptotic values of <r 2 (t)> could only 
be extrapolated in Fig. [6] The estimated diffusion 
coefficient takes the value of D a =(0.2±0.2)xI0 -5 cm 2 /s. 
for both a=Ge and Se. To allow a comparison with 
experimental data, we use the Eyring relationship; 



D 



kT 
^A 



(4) 



where A is a typical hopping length for the diffus- 
ing atomic and rj is the viscosity. In silicate melts 
such as Na,2SiC>3, Eq. Q holds satisfactorily with 
A=2.8 A, a distance typical of Si-Si and 0-0 separation 
in these melts.— By exploiting viscosity data of liquid 
GeSe2, together with a choice of A close to the Gc— Gc 
and Se— Se separation distances (3.7 A), we obtain for 
the temperatures T = 1050 K a diffusion coefficient 
-D Q ,e2;p=0.045xI0 -5 cm 2 s _1 /s4i This estimate can be 
taken as a lower bound in Eq. Q, since the account 
of homopolar bonds for both species increases A. The 
present calculations are in better agreement with the es- 
timated value of D a than the results of Ref. [H. Such 
improved agreement is consistent with the partially re- 
stored tetrahedral order that characterizes our real space 
data. Restored tetrahedral order enhances the strenght 
of the network, thereby reducing the atomic mobility. 



VI. DISCUSSION 

In view of the above pieces of evidence, it appears that 
the structure of a prototypical network-forming mate- 
rial, such as Z-GeSe2, is highly sensitive to the specific 
exchange-correlation functionals employed with DFT. In 
particular, the extent of the departure from perfect chem- 
ical order (i.e. no homopolar bonds and all Ge atoms 



four-fold coordinated) depends on the spatial distribu- 
tion of the valence charge density. In the case of moder- 
ate difference of electronegativity between the compo- 
nent systems, the account of the delicate balance be- 
tween electron localization on the atomic sites (i.e. ionic 
bonding) and electronic derealization (i.e. covalent ef- 
fects and/or tendency toward metallic bonding) becomes 
a challenging issue for atomic-scale modelling. Liquid 
GeSe2 proved to be an interesting benchmark systems for 
theoretical approaches, due to the concomitant presence 
of close percentages of corner- and edge-sharing connec- 
tions, homopolar bonds and fluctuations of concentra- 
tions on intermediate range distances. Early effective 
potentials based on the coulombic interations between 
formal charges provide the extreme case of ionic interac- 
tions, leading to networks made of undefective tetrahe- 
dra and the absence of homopolar bonds.— The inclusion 
of polarization effects and/or many body forces allows 
for the presence of edge-sharing tetrahedra, in agreement 
with experimental results ! 8 ' 9 ' 10 ? 38 The use of DFT mod- 
els proved necessary to obtain defective tetrahedra (Ge 
atoms not fourfold coordinated) and homopolar bonds 
for both Se and Ge atoms, as found in the experiments. 
However, the mere application of a fully self-consistent 
LDA had the effect of overestimating the chemical disor- 
der. This has led to the lack of intermediate range order 
and to an overall atomic structure in worse agreement 
with the experiment than that obtained with interatomic 
potentials^ These shortcomings were traced back to an 
insufficient treatment of the ionicity, i.e. an insufficient 
valence charge localization on the atomic sites, strongly 
affecting the stability of GeSe4 tetrahedra. To go be- 
yond these limitations of the LDA-DFT approach, we 
have applied the generalized gradient approximation of 
DFT in two subsequent steps. In the first (Ref. [lj) a 
functional rooted on the electron gas model has allowed 
to recover a very good agreement with the experiments in 
terms of structural properties. In the second (this work) 
a further refinement aimed at a more enhanced valence 
charge localization on the atomic site has improved the 
short range properties, in particular those related to the 
Ge close environment. 

It is useful in this context to recall the role played 
by specific structural features in determining a realis- 
tic network structure. The following ideas apply to the 
whole family of A„Xi_„ (A=Ge, Si; X=0, S, Se) disor- 
dered systems. The structural features we refer to are : 
a) homopolar bonds and miscoordinations and b) edge- 
sharing tetrahedra vs corner-sharing tetrahedra. Being 
able to correctly describe the departure from chemical or- 
der (occurring because of miscoordinations and homopo- 
lar bonds) is crucial for two reasons. First, it allows to 
model the short range structure reliably improving upon 
simplified models consisting of perfect tetrahedra con- 
nected to each other. Second, it has a profound impact 
on the intermediate range properties, which were found 
to be strongly dependent on the existence of a moderate 
amount of chemical disorder. As detailed in a previous 
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paper, no intermediate range order exist either in per- 
fect tetrahedral or in highly disordered networks^ The 
establishement of the intermediate range order is also 
very much sensitive to the occurrence of edge-sharing 
tetrahedra. Recently, it has been shown that the first 
sharp diffraction peak in the concentration-concentration 
structure factor Scc(k) is due to chains of edge-sharing 
tetrahedral This result correlates a specific structural 
subunit to the intermediate range order, exemplyifing 
how an accurate description allows to link microscopic 
features to measurable properties. 

VII. CONCLUSION 

We have studied the short and intermediate range 
properties of liquid GeSe2 at T=1050 K in the framework 
of first-principle molecular dynamics by using for the gen- 
eralized gradient approximation the BLYP scheme. Our 
motivation rests on the outcome of a previous investi- 
gation, carried out within the same theoretical frame- 
work but employing a different recipe (the one due to 
Perdew and Wang, PW) for the generalized gradient 
approximation^ In that work, it was stressed that the 
comparison between theory and experiment was excel- 
lent at the level of the total neutron structure factor, but 
less satisfactory in terms of the partial correlations both 
in real and reciprocal space. We recall that in Ref. Il5l . 
the Ge— Ge pair correlation function was less structured 



than the experimental counterpart, the first neighbors 
distances exceeded the experimental values by about 15 
% and no features were found at the FSDP location in the 
concentration-concentration structure factor. The use of 
the BLYP exchange-correlation functional substantially 
improves the short-range properties of liquid GeSe2 . The 
shapes of the Ge— Ge and (to a lesser extent) the Ge— Se 
pair correlation functional are more structured and are 
consistent with a higher level of tetrahedral organiza- 
tion, i.e. the network is less affected by coordinations 
other than the fourfold one (GeSe4 tetrahedra). The 
same occurs for the Se— Ge— Se and Ge— Se— Ge bond 
angle distributions, the first more symmetric around 
the tetrahedral angle 109° while the second shows two 
distinct peaks accounting for edge-sharing and corner- 
sharing connections. Higher tetrahedral order results in 
lower diffusion coefficient for both species. The impact 
of the BLYP scheme on the intermediate range proper- 
ties is more elusive. Smaller improvements are found for 
the intensity and the position of the peaks located at 
low k values in the partial structure factor. Work is in 
progress to investigate the impact of hybrid exchange- 
correlation functionals on the network properties of dis- 
ordered chalcogenides42 
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TABLE I: First (FPP) and second (SPP) peak positions 
in experimental (Ref. [yl3) and theoretical g a a(s)- BLYP: 
present calculations, PW: calculation of Ref. Ha . The in- 
tegration ranges corresponding to the coordination numbers 
n a j3 and n' a p are 0—2.6 A, 2.6—4.2 A for gGeGe(r), 0—3.1 A, 
3.1-4.5 A for g G eSe{r) and 0-2.7 A, 2.7-4.8 A for g S eSe{i)- 
Error bars are the standard deviations from the mean for sub- 
averages of 2 ps. 



TABLE II: Experimental and theoretical values for the partial 
coordination numbers nGe and nse and the average coordina- 
tion number n of liquid GeSe2 at T=1040 K. BLYP: present 
calculations, PW: calculation of Ref. [l5|. The coordination 
numbers n Ge and nse are given by nceGe + nceSe and nseSe 
+ nseGe, respectively (see the values reported in Table U for 
riGeGe, riGeSe and nseSe, where riGeSe = 2ns eGe)- The aver- 
age coordination number n is equal to CGe(n,GeGe + ^GeSe) 
+ cs e (nse.Se + fiSeGe). The experimental values extracted 
from Ref. are also reported. Error bars are the standard 
deviations of the mean for subaverages of 2 ps. 



BLYP 
PW 
Ref. 6 



U G e 

3.77±0.02 
3.80±0.02 
3.75±0.3 



nse 

2.11±0.02 
2.25±0.02 
1.98±0.15 



n 

2.66±0.02 
2.77±0.02 
2.57±0.20 



exp / \ 
9 G eGe( r ) 

Sg57(r) 
<&(r) 

9GeSe( T ) 

sfX(r) 

exp i \ 
dSeSei 1 ) 



FPP (A) 

2.45±0.10 
2.70±0.10 
2.33±0.03 

2.36±0.10 
2.41±0.10 
2.42±0.02 

2.38±0.02 
2.34±0.02 
2.30±0.02 



n a p 

0.22±0.01 
0.04±0.01 
0.25±0.10 

3.55±0.01 
3.76±0.01 
3.5±0.2 

0.33±0.01 
0.37±0.01 
0.23±0.05 



SPP (A) 

3.67±0.10 
3.74±0.05 
3.59±0.02 

5.67±0.02 
5.60±0.01 
4.15±0.10 

3.83±0.02 
3.84±0.02 
3.75±0.02 



n a {3 

2.70±0.06 
2.74±0.06 
2.9±0.3 

3.85±0.06 
3.72±0.03 
4.0±0.3 

8.9±0.06 

9.28±0.04 

9.6±0.3 
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TABLE III: Average number n a (l) (expressed as a percent- 
age) of Ge and Se atoms Z-fold coordinated at a distance of 3.0 
A. For each value of n a (l), we give the identity and the num- 
ber of the Ge and Se neighbors. For instance, GeSe3 with ?=4 
means a fourfold coordinated Ge with one Ge and three Se 
nearest-neighbors. Values smaller than 1 are reported only for 
sake of comparison with corresponding values equal or larger 
than 1. In parenthesis: results of Ref. ITsl . We also com- 
pare calculated and experimental values (in percentage) for 
the number of Ge atoms forming edge-sharing connections, 
./Vgc(ES), the number of Ge atoms forming corner-sharing 
connections, iVGc(CS) and the number of Ge atoms involved 
in homopolar bonds, A^gb-Gg- Experimental values are taken 
from Ref. & 
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I = 2 
Se 2 


4.0 (5.2) 


I = 3 

GeSe 2 

Se 3 


0.8 (2.6) 
13.5 (19.8) 




I = 4 

GeSe 3 

Se 4 


23.3 (7.0) 
41.8 (53.8) 


I = 5 
Ge 2 Se3 
GeSe 4 
Se 5 


2.4 (0.4) 
11.7 (5.9) 
0.6 (4.6) 


Se 


I = 1 
Ge 


0.9 (1.7) 


I = 2 
Se 2 
SeGe 
Ge 2 


3.6 (2.8) 
20.7 (21.9) 
59.2 (45.6) 




I = 3 
Se 2 Ge 
SeGe 2 
Ge 3 


2.4 (3.2) 
5.8 (8.6) 
6.4 (13.1) 


I = 4 
SeGe 3 
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This work 
Ref. 15 


N Ge {ES) 

39 

45 


iVGc(CS) 

39 

51 


A^Ge-Ge 

22 
4 





Energy (eV) 



FIG. 1: (color on line) Electronic density of states (Kohn- 
Sham eigenvalues) of liquid GeSe2i black line: present BLYP 
results, red line: PW results of Ref. [l5|. A gaussian broaden- 
ing of 0.1 eV has been employed. 




FIG. 2: (color on line) Partial pair correlation functions for 
liquid GeSe2: black line: present BLYP results, red line: PW 
results of Ref. [T5 , open circles: experimental results of Ref. 0- 




FIG. 3: (color on line) Bond-angle distributions Ge— Se— Ge 
(top) and Se— Ge— Se (bottom), black line: present BLYP 
results, red line: PW results of Ref. [la . 
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FIG. 4: Faber-Ziman partial structure factors for liquid 
GeSe2- In each panel, the experimental results of Ref. are 
compared with the present BLYP results (bottom part) and 
with the P W results of Ref. Gl (top part) . Sa% e W , Sglse (&) 
and Ss^s e (k) have been shifted up by 4, 2, and 2, respectively. 
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k(A"') 

FIG. 5: Bhatia- Thornton partial structure factors for liquid 
GeSe2- In each panel, the experimental results of Ref. are 
compared with the present BLYP calculations (bottom part) 
and with the PW calculations of Ref. Gl (top part). Sf?g(k), 
Snc (k) an d Sec (k) have been shifted up by 1, 0.6, and 0.5, 
respectively. 
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FIG. 6: Average mean square displacements for Ge (dashed 
line) and Se atoms (full line) of liquid GeSe2. The infinite time 
behavior corresponding to a slope equal to 1 in the log— log 
plot of <r 2 (t)> vs t is also shown. 



